
function [c,ceq,Gc,Gceq] = Shares(B,S,X,KD,A,n,nu,wnu)

N0 = n(1); N1 = n(2); N2 = n(3); N = N0 + N1 + N2;
[NN,K] = size(X);
KK = length(B) - (K + 2 + NN);
R = size(nu,1);

c = []; % inequality constraints

s = X * B(1:K,1) + B(K+1,1) * X(:,K-1) * nu(:,1)' + ...
    B(K+2,1) * X(:,K) * nu(:,2)' + B(K+2+(1:NN),1);
s = exp(s);
s = s ./ (1 + A * s); % predicted vote shares

ceq = s * wnu - S; % equality constraints

if nargout > 2
    
    Gc = []; % Jacobian of inequality constraints
    Gceq = zeros(NN,K+2+NN+KK); % Jacobian of equality constraints
    ss = zeros(N,R,5);
    ss(:,:,1) = s(1:N,:);
    ss(:,:,2) = s(N+(1:N),:);
    ss(1:N0,:,3) = s(2*N+(1:N0),:);
    ss(N0+N1+(1:N2),:,3) = s(2*N+N0+(1:N2),:);
    ss(1:N0+N1,:,4) = s(2*N+N0+N2+(1:N0+N1),:);
    ss(:,:,5) = s(3*N+N0+(1:N),:);
    xx = zeros(N,1,5);
    for k = 1:K
        xx(:,:,1) = X(1:N,k);
        xx(:,:,2) = X(N+(1:N),k);
        xx(1:N0,:,3) = X(2*N+(1:N0),k);
        xx(N0+N1+(1:N2),:,3) = X(2*N+N0+(1:N2),k);
        xx(1:N0+N1,:,4) = X(2*N+N0+N2+(1:N0+N1),k);
        xx(:,:,5) = X(3*N+N0+(1:N),k);
        sx = sum(ss.*xx,3);
        Gceq(1:N,k) = (ss(:,:,1) .* (xx(:,:,1) - sx)) * wnu;
        Gceq(N+(1:N),k) = (ss(:,:,2) .* (xx(:,:,2) - sx)) * wnu;
        Gceq(2*N+(1:(N0+N2)),k) = (ss([1:N0,N0+N1+(1:N2)],:,3) .* (xx([1:N0,N0+N1+(1:N2)],:,3) ...
            - sum(ss([1:N0,N0+N1+(1:N2)],:,:) .* xx([1:N0,N0+N1+(1:N2)],:,:),3))) * wnu;
        Gceq(2*N+N0+N2+(1:N0+N1),k)=(ss(1:N0+N1,:,4) .* (xx(1:N0+N1,:,4) ...
            - sum(ss(1:N0+N1,:,:) .* xx(1:N0+N1,:,:),3))) * wnu;
        Gceq(3*N+N0+(1:N),k) = (ss(:,:,5) .* (xx(:,:,5) - sx)) * wnu;
    end
    xx = zeros(N,R,5);
    for k = 1:2
        xx(:,:,1) = X(1:N,K-2+k) * nu(:,k)';
        xx(:,:,2) = X(N+(1:N),K-2+k) * nu(:,k)';
        xx(1:N0,:,3) = X(2*N+(1:N0),K-2+k) * nu(:,k)';
        xx(N0+N1+(1:N2),:,3) = X(2*N+N0+(1:N2),K-2+k) * nu(:,k)';
        xx(1:N0+N1,:,4) = X(2*N+N0+N2+(1:N0+N1),K-2+k) * nu(:,k)';
        xx(:,:,5) = X(3*N+N0+(1:N),K-2+k) * nu(:,k)';
        sx = sum(ss.*xx,3);
        Gceq(1:N,K+k) = (ss(:,:,1) .* (xx(:,:,1) - sx)) * wnu;
        Gceq(N+(1:N),K+k) = (ss(:,:,2) .* (xx(:,:,2) - sx)) * wnu;
        Gceq(2*N+(1:(N0+N2)),K+k) = (ss([1:N0,N0+N1+(1:N2)],:,3) .* (xx([1:N0,N0+N1+(1:N2)],:,3) ...
            - sum(ss([1:N0,N0+N1+(1:N2)],:,:) .* xx([1:N0,N0+N1+(1:N2)],:,:),3))) * wnu;
        Gceq(2*N+N0+N2+(1:N0+N1),K+k) = (ss(1:N0+N1,:,4) .* (xx(1:N0+N1,:,4) ...
            - sum(ss(1:N0+N1,:,:) .* xx(1:N0+N1,:,:),3))) * wnu;
        Gceq(3*N+N0+(1:N),K+k) = (ss(:,:,5) .* (xx(:,:,5) - sx)) * wnu;
    end
    for t = 1:N0
        tt = [t,N+t,2*N+t,2*N+N0+N2+t,3*N+N0+t];
        Gceq(tt,K+2+tt) = Gceq(tt,1+[0,KD+3,2*(KD+3),2*(KD+3)+KD+2,2*(KD+3)+2*(KD+2)]);
    end
    for t = 1:N1
        tt = [N0+t,N+N0+t,2*N+N0+N2+N0+t,3*N+N0+N0+t];
        Gceq(tt,K+2+tt) = Gceq(tt,2+[0,KD+3,2*(KD+3)+KD+2,2*(KD+3)+2*(KD+2)]);
    end
    for t = 1:N2
        tt = [N0+N1+t,N+N0+N1+t,2*N+N0+t,3*N+N0+N0+N1+t];
        Gceq(tt,K+2+tt) = Gceq(tt,3+[0,KD+3,2*(KD+3)-1,2*(KD+3)+2*(KD+2)]);
    end
    
    Gceq = sparse(Gceq');

end



